Statistical analyses and plotting for sorghum rhizosphere GWAS
Siwen Deng Ph.D. and Daniel F. Caddell Ph.D.
2/24/2020
all R code for the manuscript entitled ‘Genome wide association study reveals plant loci controlling heritability of the rhizosphere microbiome’
Assume the input files below are in the current working directory
load packages
library("phyloseq"); packageVersion("phyloseq") #‘1.30.0’
library("ggplot2"); packageVersion("ggplot2") #‘3.2.1’
library("scales"); packageVersion("scales") #‘1.1.0’
library("grid"); packageVersion("grid") #‘3.6.1’
library("DESeq"); packageVersion("DESeq") #‘1.38.0’
library("ape"); packageVersion("ape") #‘5.3’
library("reshape2"); packageVersion("reshape2") #‘1.4.3’
library("vegan"); packageVersion("vegan") #‘2.5.6’
library("data.table"); packageVersion("data.table") #‘1.12.8’
library("ggrepel"); packageVersion("ggrepel") #‘0.8.1’
library("dplyr"); packageVersion("dplyr") #‘0.8.3’
library("data.table"); packageVersion("data.table") #‘1.12.8’
library("qqman"); packageVersion("qqman") #‘0.1.4’
library("RColorBrewer"); packageVersion("RColorBrewer") #‘1.1.2’
library("colorspace"); packageVersion("colorspace") #‘1.4.1’
library("patchwork"); packageVersion("patchwork") #‘1.0.0’
theme_set(theme_bw())
#setwd("/set/working/directory/location/")
setwd("/users/dcaddell/Desktop/Github/R_data/")Figure 1
Figure 1b and figure 1c
#import data from current working directory
rar <- readRDS("fig1_24line.rds")
Diversity_table <- estimate_richness(rar, measures=c("Shannon", "Observed"))
Diversity_table <- merge(Diversity_table, sample_data(rar), by="row.names")
Diversity_table$Sampletype <- factor(Diversity_table$Sampletype, levels = c( "Leaf", "Root", "Rhizo"))Plot Figure 1b and Figure 1c
p1b <- ggplot(data=Diversity_table, aes(y=Shannon,x=Sampletype,fill=Sampletype)) +
geom_boxplot() +
scale_fill_manual(values=c("#9ACAA1","#E1D337", "#CE8764")) +
theme(axis.text.x=element_text(hjust=1,vjust=0.5,size=10,color="black",angle=90,face="bold"),
axis.text.y=element_text(size=11,color="black",face="bold"),
axis.title=element_text(size=11,face="bold"),text=element_text(size=11,face="bold"),
legend.position="right")
p1b
# ggsave("figure_1b.pdf", plot = p1b, width=6, height=5, useDingbats=FALSE)
# ggsave("figure_1b.png", plot = p1b, width=6, height=5, dpi=600)
p1c <- plot_ordination(rar, ordinate(rar, "MDS",distance="bray"),axes=1:2, color = "Sampletype") +
scale_color_manual(values=c("#9ACAA1","#CE8764","#E1D337")) +
geom_point(colour="black",size=3.5) +
geom_point(size = 2.5) +
xlim(-0.4,0.7) +
scale_x_continuous(breaks=c(-0.3,0,0.3,0.6))+
ylim(-0.33,0.25)+
theme(axis.text.x=element_text(size=11,color="black",angle=90),
axis.text.y=element_text(size=11,color="black"),
axis.title=element_text(size=11,face="bold"),
text=element_text(size=11,face="bold"),
legend.position="right")## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.p1c
# ggsave("figure_1c.pdf", plot = p1c, width=6, height=5, useDingbats=FALSE)
# ggsave("figure_1c.png", plot = p1c, width=6, height=5, dpi=600)Figure 1d
# to create Root subset
root <- subset_samples(rar, Sampletype == "Root")
raw_env <- data.frame(sample_data(root))
raw_env$sbNum_rep<-paste(raw_env$Line,raw_env$Block,sep="_")
sample_data(root) <- raw_env
OTU1 <- data.frame(otu_table(root))
rownames(OTU1) <- sub("^", "X", rownames(OTU1))
OTU1 = as(OTU1, "matrix")
OTU1 <- t(OTU1)
Root_OTU = as.data.frame(OTU1)
env_current <- data.frame(sample_data(root))
row.names(env_current)<-env_current$sbNum_rep
row.names(Root_OTU)<-env_current$sbNum_rep
Root_OTU<-Root_OTU[order(row.names(Root_OTU)),]
# to create Rhizo subset
rhizo <- subset_samples(rar, Sampletype == "Rhizo")
raw_env <- data.frame(sample_data(rhizo))
raw_env$sbNum_rep<-paste(raw_env$Line,raw_env$Block,sep="_")
sample_data(rhizo) <- raw_env
OTU1 <- data.frame(otu_table(rhizo))
rownames(OTU1) <- sub("^", "X", rownames(OTU1))
OTU1 = as(OTU1, "matrix")
# transpose if necessary
OTU1 <- t(OTU1)
# Coerce to data.frame
Rhizo_OTU = as.data.frame(OTU1)
env_current <- data.frame(sample_data(rhizo))
row.names(env_current)<-env_current$sbNum_rep
row.names(Rhizo_OTU)<-env_current$sbNum_rep
Rhizo_OTU<-Rhizo_OTU[order(row.names(Rhizo_OTU)),]
# to create the Leaf subset
leaf <- subset_samples(rar, Sampletype == "Leaf")
raw_env <- data.frame(sample_data(leaf))
raw_env$sbNum_rep<-paste(raw_env$Line,raw_env$Block,sep="_")
sample_data(leaf) <- raw_env
OTU1 <- data.frame(otu_table(leaf))
rownames(OTU1) <- sub("^", "X", rownames(OTU1))
OTU1 = as(OTU1, "matrix")
# transpose if necessary
OTU1 <- t(OTU1)
# Coerce to data.frame
Leaf_OTU = as.data.frame(OTU1)
env_current <- data.frame(sample_data(leaf))
row.names(env_current)<-env_current$sbNum_rep
row.names(Leaf_OTU)<-env_current$sbNum_rep
Leaf_OTU<-Leaf_OTU[order(row.names(Leaf_OTU)),]
#to create the distance objects for each subset
Root_dist<-vegdist(Root_OTU,"bray")
Rhizo_dist<-vegdist(Rhizo_OTU,"bray")
Leaf_dist<-vegdist(Leaf_OTU,"bray")
## to create the phylogenetic distance table for each OTU subset
phy_distances_table<-read.table("fig1_24line.txt")
phy_distances_table_Root<-phy_distances_table[c(row.names(Root_OTU)),c(row.names(Root_OTU))]
phy_distances_table_Rhizo<-phy_distances_table[c(row.names(Rhizo_OTU)),c(row.names(Rhizo_OTU))]
phy_distances_table_Leaf<-phy_distances_table[c(row.names(Leaf_OTU)),c(row.names(Leaf_OTU))]
#to create the phylogenetic distance for each distance table
phy_dist_Root<-as.dist(phy_distances_table_Root)
phy_dist_Rhizo<-as.dist(phy_distances_table_Rhizo)
phy_dist_Leaf<-as.dist(phy_distances_table_Leaf)Perform Mantels test with vegan and stores values in a new dataframe
test<-character()
Mstat<-numeric()
signif<-numeric()
mantels<-data.frame(test,Mstat,signif)
#spearman
mantel_Root<-vegan::mantel(phy_dist_Root, Root_dist, method="spearman",permutations=9999,strata=NULL)
mantels<-rbind(mantels, data.frame(test=c("phy_Root"),Mstat=mantel_Root$statistic,signif=mantel_Root$signif))
mantel_Rhizo<-vegan::mantel(phy_dist_Rhizo, Rhizo_dist, method="spearman",permutations=9999,strata=NULL)
mantels<-rbind(mantels, data.frame(test=c("phy_Rhizo"),Mstat=mantel_Rhizo$statistic,signif=mantel_Rhizo$signif))
mantel_Leaf<-vegan::mantel(phy_dist_Leaf, Leaf_dist, method="spearman",permutations=9999,strata=NULL)
mantels<-rbind(mantels, data.frame(test=c("phy_Leaf"),Mstat=mantel_Leaf$statistic,signif=mantel_Leaf$signif))
mantels<-cbind(mantels,reshape2::colsplit(mantels$test,pattern="_",names=c("None","SampleType")))
mantels$SampleType<-factor(mantels$SampleType, levels=c("Leaf","Root","Rhizo"))
mantels## test Mstat signif None SampleType
## 1 phy_Root 0.06136486 0.1223 phy Root
## 2 phy_Rhizo 0.12999624 0.0178 phy Rhizo
## 3 phy_Leaf -0.03982383 0.6827 phy Leaf# write.table(mantels, "figure1d_mantel.csv", sep=",", row.names=FALSE, quote=FALSE)Plot Figure 1d
p1d <- ggplot(data=mantels,aes(x=SampleType,y=Mstat, fill=SampleType)) +
geom_bar(stat="identity")+
scale_fill_manual(values=c("#9ACAA1","#E1D337", "#CE8764"))+
ylim(-0.05,0.15)
p1d
# ggsave("figure_1d.pdf", plot = p1d, width=6, height=5, useDingbats=FALSE)
# ggsave("figure_1d.png", plot = p1d, width=6, height=5, dpi=600)Figure 2
#Figure 2 color palette
F2B_palette <- c( "#599ec4" #ST Giles Blue
,"#a1c5c8" #Blue Ground
,"#7997a1" #Stone Blue
,"#427e83" #Vardo
,"#84b59c" #Arsenic
,"#919f70" #Yeabridge Green
,"#686a47" #Bancha
,"#c8bd83" #Churlish Green
,"#cb9e59" #India Yellow
,"#ecc363" #Babouche
,"#c57b67" #Red Earth
,"#d65f3d" #Charlotte's Locks
,"#a04344" #Incarnadine
,"#bf7a8f" #Rangwali
,"#8B0000"#"dark red"
,"#8d8089" #Brassica
,"#e5e0db" #Strong White
,"#444546") #Off-Black
#Top 17 bacterial orders sorted by heritability (at least 6 OTUs in an order)
sort_order <- c("Sphingobacteriales",
"Actinomycetales",
"Verrucomicrobiales",
"Myxococcales",
"Gemmatimonadales",
"Rhodospirillales",
"Planctomycetales",
"Burkholderiales",
"Acidobacteriales",
"Rhizobiales",
"Solirubrobacterales",
"Xanthomonadales",
"Solibacterales",
"Spartobacteriales",
"Bacillales",
"TM7PH",
"Flavobacteriales",
"Other")Figure 2a
#import data from current working directory
df <- read.table("fig2_count.txt", header = TRUE)
df$Order <- factor(df$Order, levels = rev(sort_order))
df$Group <- factor(df$Group, levels = c("all", "non_herit", "herit"))Plot Figure 2a
p2a <- ggplot(data=df, aes(x=Group, y=Count, fill=factor(Order))) +
geom_bar(stat="identity",position = "fill") +
scale_fill_manual(values = rev(F2B_palette)) +
theme_minimal()+
theme(axis.text.x=element_text(size=12,color="black",angle=270,hjust=0, vjust=0.5),
axis.text.y=element_text(size=12,color="black"),
axis.title=element_text(size=16,face="bold"),
text=element_text(size=16))
p2a
# ggsave("figure_2a.pdf", plot = p2a, width=4.5, height=6, useDingbats=FALSE)
# ggsave("figure_2a.png", plot = p2a, width=4.5, height=6, dpi=600)Figure 2b
#import data from current working directory
df <- read.table("fig2_relative.txt", header = TRUE)
df$Top17_color <- factor(df$Top17_color, levels = sort_order)Plot Figure 2b
p2b <- ggplot() +
geom_point(data=df, mapping=aes(x=NormCount, y=NormRelAbun, color=Top17_color, size=all_abun)) +
theme_minimal() +
xlim(-4,4) +
ylim(-5,5)+
labs(x = "OTU counts \n Log2(heritable/non-heritable)", y = "Read count abundance \n Log2(heritable/non-heritable)") +
scale_color_manual(values = F2B_palette) +
geom_vline(xintercept = 0) +
geom_hline(yintercept = 0)
p2b
# ggsave("figure_2b.pdf", plot = p2b, width=8, height=6, useDingbats=FALSE)
# ggsave("figure_2b.png", plot = p2b, width=8, height=6, dpi=600)Figure 3
#subset overlapping orders
S_M10_M15_only15 <- c("Caulobacterales","Burkholderiales","Sphingobacteriales","Actinomycetales","Flavobacteriales","Sphingomonadales","Myxococcales","Bacillales","Xanthomonadales","Gemmatimonadales","Verrucomicrobiales","Rhodospirillales","Rhizobiales","Pseudomonadales","BetaproteobacteriaCL")
m10_m15_all20 <- c("Caulobacterales","Burkholderiales","Sphingobacteriales","Actinomycetales","Flavobacteriales","Sphingomonadales","Myxococcales","Bacillales","Xanthomonadales","Gemmatimonadales","Verrucomicrobiales","Rhodospirillales","Rhizobiales","Pseudomonadales","BetaproteobacteriaCL","Chitinophagales","Pedosphaerales","Nitrospirales","Nevskiales","MND1")
m10_only17 <- c("Caulobacterales","Burkholderiales","Sphingobacteriales","Actinomycetales","Flavobacteriales","Sphingomonadales","Myxococcales","Bacillales","Xanthomonadales","Gemmatimonadales","Verrucomicrobiales","Rhodospirillales","Rhizobiales","Pseudomonadales","BetaproteobacteriaCL","Solibacterales","A4b")
m10_all28 <- c("Caulobacterales","Burkholderiales","Sphingobacteriales","Actinomycetales","Flavobacteriales","Sphingomonadales","Myxococcales","Bacillales","Xanthomonadales","Gemmatimonadales","Verrucomicrobiales","Rhodospirillales","Rhizobiales","Pseudomonadales","BetaproteobacteriaCL","Solibacterales","A4b","Nitrososphaerales","SC-I-84","258ds10","Opitutales","AcidobacteriaCL","Micromonosporales","Chitinophagales","Pedosphaerales","Nitrospirales","Nevskiales","MND1")
m15_only19 <- c("Caulobacterales","Burkholderiales","Sphingobacteriales","Actinomycetales","Flavobacteriales","Sphingomonadales","Myxococcales","Bacillales","Xanthomonadales","Gemmatimonadales","Verrucomicrobiales","Rhodospirillales","Rhizobiales","Pseudomonadales","BetaproteobacteriaCL","Planctomycetales","Acidobacteriales","Methylophilales","Thiotrichales")
m15_all47 <- c("Xanthomonadales", "Gemmatimonadales", "Verrucomicrobiales", "Rhodospirillales", "Rhizobiales", "Pseudomonadales", "BetaproteobacteriaCL", "Planctomycetales", "Acidobacteriales", "Methylophilales", "Thiotrichales", "Gemm-3CL", "MIZ46", "Solirubrobacterales", "Ellin6529CL", "Sva0725", "Gaiellales", "Acidimicrobiales", "Ellin5290", "Bdellovibrionales", "Rhodobacterales", "Gemm-1CL", "Rickettsiales", "AKYG885", "C0119CL", "Pirellulales", "Bacteroidales", "order 11-24", "0319-7L14", "Spirobacillales", "Micrococcales", "Legionellales", "Unassigned", "Ellin6067", "Chitinophagales", "Pedosphaerales", "Nitrospirales", "Nevskiales", "MND1")
#import data from current working directory
rar1 <- readRDS("fig3_200line.rds")
otus <- taxa_names(rar1)
out <- data.frame(subset=character(),
S_M10_M15_only15=numeric(),
m10_m15_all20=numeric(),
m10_only17=numeric(),
m10_all28=numeric(),
m15_only19=numeric(),
m15_all47=numeric(),
stringsAsFactors=FALSE)
#perform subsampling
for(i in 1:10000){
print(i)
subset <- sample(otus, 100, replace = FALSE, prob = NULL)
all_taxa <- as.data.frame(rar1@tax_table@.Data)
subset_taxa <- subset(all_taxa, row.names(all_taxa) %in% subset)
subset_order <- unique(factor(subset_taxa$Order))
a <- length(intersect(subset_order, S_M10_M15_only15))
b <- length(intersect(subset_order, m10_m15_all20))
c <- length(intersect(subset_order, m10_only17))
d <- length(intersect(subset_order, m10_all28))
e <- length(intersect(subset_order, m15_only19))
f <- length(intersect(subset_order, m15_all47))
add <- data.frame(subset=paste0("subset",1),
S_M10_M15_only15=a,
m10_m15_all20=b,
m10_only17=c,
m10_all28=d,
m15_only19=e,
m15_all47=f)
out <- rbind(out, add)
}Export results of 10,000 permutations
plot(density(out$S_M10_M15_only15))
plot(density(out$m10_m15_all20))
plot(density(out$m10_only17))
plot(density(out$m10_all28))
plot(density(out$m15_only19))
plot(density(out$m15_all47))
# write.table(out, "figure_3_permutations.csv", sep=",", row.names=FALSE, quote=FALSE)Figure 4
#import data from current working directory
qval <- read.csv("fig4_pc1.txt", header = TRUE, sep = "\t")
qval$log10p <- -log10(qval$p_score)
geno <- data.frame(CHR=qval$chr, BP=qval$ps, P=qval$p_score, SNP=qval$rs)Plot figure 4a
p4a <- qqman::manhattan(na.omit(geno), suggestiveline = F, col = c("#E09E19", "#2D637F")) #export pdf size width=10, height=3.5, png width = 1000, height = 350
Figure 4b allele ratio heatmap
#import data from current working directory
ch4_figure4 <- read.csv("fig4_ch4.txt")
#Figure 4b color palette
F4B_palette <- c("#cb9e59"#"Acidobacteriales"
,"#a1c5c8"#"Actinomycetales"
,"#8B0000"#"Bacillales"
,"#c8bd83"#"Burkholderiales"
,"#84b59c"#"Gemmatimonadales"
,"#427e83"#"Myxococcales"
,"#444546" #"Other"
,"#a04344"#"Solibacterales"
,"#c57b67"#"Solirubrobacterales"
,"#bf7a8f"#"Spartobacteriales"
,"#599ec4"#"Sphingobacteriales"
,"#7997a1"#"Verrucomicrobiales"
,"#d65f3d")#"Xanthomonadales" Plot figure 4b allele ratio heatmap
#plot major/minor allele ratio (Log2) with red/blue gradient
p4b1 <- ggplot(data=ch4_figure4, aes(x=Manual_Rank_fig4, y=1, fill= Log2ratio)) +
geom_tile(color = "white")+
scale_fill_gradient2(low = "#4f94cd", high = "#ee2c2c", mid = "white",
midpoint = 0, limit = c(-1.5,1), space = "Lab") +
theme_minimal()+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, size = 12, hjust = 1),
axis.title.y=element_blank())+
coord_fixed()
p4b1
# ggsave("figure_4b-1.pdf", plot = p4b1, width=10, height=3, useDingbats=FALSE)
# ggsave("figure_4b-1.png", plot = p4b1, width=10, height=3, dpi=600)
#plot monoderm/diderm
p4b2 <- ggplot(data=ch4_figure4, aes(x=Manual_Rank_fig4, y=1, fill= MonoDiderm)) +
geom_tile()+
theme_minimal()+
scale_fill_manual(values = c("#50414c","#e5e0db","#cb9e59" )) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, size = 12, hjust = 1),
axis.title.y=element_blank())+
coord_fixed()
p4b2
# ggsave("figure_4b-2.pdf", plot = p4b2, width=10, height=3, useDingbats=FALSE)
# ggsave("figure_4b-2.png", plot = p4b2, width=10, height=3, dpi=600)
#plot orders
p4b3 <- ggplot(data=ch4_figure4, aes(x=Manual_Rank_fig4, y=1, fill= Order_fig4)) +
geom_tile()+
theme_minimal()+
scale_fill_manual(values = (F4B_palette)) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, size = 12, hjust = 1),
axis.title.y=element_blank())+
coord_fixed()
p4b3
# ggsave("figure_4b-3.pdf", plot = p4b3, width=10, height=6, useDingbats=FALSE)
# ggsave("figure_4b-3.png", plot = p4b3, width=10, height=6, dpi=600)Figure 4b GWAS heatmap
#import data from current working directory
goi <- readRDS("fig4_ch4.rds")
##rank OTUs
levels_otu <- c(rev(c("X1039","X172","X544","X798","X1170","X78","X30","X506","X996",
"X1289","X11","X325","X1239","X326","X461","X904","X943","X190",
"X1020","X199","X1190","X302","X692","X14","X695","X286","X94","X82",
"X766","X312","X531","X734","X205","X148","X340","X318","X577","X18",
"X947","X445")),"PC1")
goi$otu <- factor(x = goi$otu, levels = levels_otu, ordered = TRUE)
length(levels(factor(goi$rs)))## [1] 64Plot figure 4b GWAS heatmap
p4b4 <- ggplot(data = goi, aes(x = rs, y = otu)) +
geom_tile(aes(fill = -log10(goi$p_score))) +
scale_fill_gradientn(colours=c("white","#e5e0db","#6a90b4","#2D637F","#4d5b6a"), values = rescale(c(0,2,4,6,8)),guide = "colorbar") +
theme(axis.title.y = element_blank(),
axis.ticks.y = element_blank(),
legend.position = "top")
p4b4
# ggsave("figure_4b-4.pdf", plot = p4b4, width=6, height=6, useDingbats=FALSE)
# ggsave("figure_4b-4.png", plot = p4b4, width=6, height=6, dpi=600)Figure 4c
#import data from current working directory
SorghumGenes<-read.table(file = "fig4_sorghum.txt",header = T,sep = "\t")
SorghumGenes <-reshape2::melt(SorghumGenes)## Using Gene_annotation as id variablesSorghumGenes = data.frame(SorghumGenes)
#rank sorghum genes in chromosome 4 locus
Gene_Anno <- c("Uncharacterized protein1",
"Disease resistance protein (RGA2)",
"chromatin assembly factor 1 subunit A (CHAF1A)",
"stearoyl-[acyl-carrier-protein] 9-desaturase 5 (SAD5)",
"F-box/kelch-repeat protein",
"exocyst complex component (EXO70B1)",
"Uncharacterized protein2",
"Uncharacterized protein3",
"beta-sesquiphellandrene synthase-like protein",
"Endo-1,4-beta-xylanase like protein",
"Uncharacterized protein4",
"Ser/Thr-protein phosphatase 6 regulatory ankyrin repeat subunit C (ANKRD52)",
"ubiquitin-like modifier-activating enzyme 5(UBA5)",
"chaperone protein, DnaJ-like",
"protein indeterminate-domain 16 (IDD16)",
"choline/ethanolamine kinase",
"Uncharacterized protein5",
"Uncharacterized protein6",
"NDR1/HIN1-like protein 6 (NHL6)",
"gamma carbonic anhydrase-like 2 (GCAL2)",
"Uncharacterized protein7",
"Uncharacterized protein8",
"Uncharacterized protein9",
"DnaJ subfamily B member",
"U6 snRNA-associated Sm-like protein (LSM6)",
"acetolactate synthase 1 (ALS1)",
"heavy metal-associated isoprenylated plant protein 7 (HIP7)")Plot figure 4c
p4c <- ggplot(data = SorghumGenes, aes(x = rev(variable), y = Gene_annotation)) +
geom_tile(aes(fill = value)) +
scale_y_discrete(limits = rev(Gene_Anno))+
scale_fill_gradientn(colours=c("white","#e5e0db","#6a90b4","#2D637F","#4d5b6a"), values = rescale(c(0,1,2,10,30)),guide = "colorbar") +
theme(axis.title.y = element_blank(),
axis.ticks = element_blank(),
axis.text.x = element_text(size=5),
axis.text.y = element_text(size=5),
legend.text=element_text(size=5))
p4c
# ggsave("figure_4c.pdf", plot = p4c, width=5, height=6, useDingbats=FALSE)
# ggsave("figure_4c.png", plot = p4c, width=5, height=6, dpi=600)Figure 5
#import data from current working directory
rar1 <- readRDS("fig5_validation.rds")
rar1## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 906 taxa and 67 samples ]
## sample_data() Sample Data: [ 67 samples by 16 sample variables ]
## tax_table() Taxonomy Table: [ 906 taxa by 7 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 906 tips and 905 internal nodes ]OTU1 = data.frame(otu_table(rar1))
rownames(OTU1) <- sub("^", "X", rownames(OTU1))
OTU1 = as(OTU1, "matrix")
# transpose if necessary
#if(taxa_are_rows(current)){OTU1 <- t(OTU1)}
# Coerce to data.frame
OTUdf <- as.data.frame(OTU1)
## orgnize metadata file
metacurrent <- data.frame(sample_data(rar1))
#rhizosphere subset
rhizo <- subset_samples(rar1, Sampletype == "Rhizosphere")
rhizo <- prune_taxa(taxa_sums(rhizo) > 0, rhizo)
rhizo## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 906 taxa and 36 samples ]
## sample_data() Sample Data: [ 36 samples by 16 sample variables ]
## tax_table() Taxonomy Table: [ 906 taxa by 7 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 906 tips and 905 internal nodes ]#perform CAP ordination
ordcap = ordinate(rhizo, "CAP", "bray", ~Group)Plot figure 5a
p5a <- plot_ordination(rhizo, ordcap, "samples", color="Group") +
geom_point(size=5) +
scale_color_manual(values=c("#a04344", "#2D637F")) +
theme(legend.title = element_text(colour="black", size=11, face="bold")) +
theme(legend.text = element_text(colour="black", size = 11, face = "bold"))+
theme(axis.text.x=element_text(hjust=1,vjust=0.5,size=10,color="black",angle=90,face="bold"),
axis.text.y=element_text(size=11,color="black",face="bold"),
axis.title=element_text(size=11,face="bold"),text=element_text(size=11,face="bold"),
legend.position = "right")
p5a
# ggsave("figure_5a.pdf", plot = p5a, width=5, height=5, useDingbats=FALSE)
# ggsave("figure_5a.png", plot = p5a, width=5, height=5, dpi=600)Figure 5b
#import data from current working directory
indic <- read.table("fig5_indicators.txt", header = TRUE)
#Figure 5b color palette
F5B_palette <-c("#cb9e59" #"Acidobacteriales"
,"#a1c5c8" #"Actinomycetales"
,"#c8bd83" #"Burkholderiales"
,"#84b59c" #"Gemmatimonadales"
,"#444546" #"Other"
,"#919f70" #"Rhodospirillales"
,"#c57b67" #"Solirubrobacterales"
,"#599ec4" #"Sphingobacteriales"
,"#d65f3d") #"Xanthomonadales"Plot figure 5b
#plot orders
p5b1 <- ggplot(data=indic, aes(x=Manual_Rank_fig5, y=1, fill= Order_fig5)) +
geom_tile()+
theme_minimal()+
scale_fill_manual(values = (F5B_palette)) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, size = 12, hjust = 1),
axis.title.y=element_blank())+
coord_fixed()
p5b1
# ggsave("figure_5b-1.pdf", plot = p5b1, width=10, height=4, useDingbats=FALSE)
# ggsave("figure_5b-1.png", plot = p5b1, width=10, height=4, dpi=600)
#plot monoderm/diderm
p5b2 <- ggplot(data=indic, aes(x=Manual_Rank, y=1, fill= CellWall)) +
geom_tile()+
theme_minimal()+
scale_fill_manual(values = c("#50414c","#e5e0db","#cb9e59" )) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, size = 12, hjust = 1),
axis.title.y=element_blank())+
coord_fixed()
p5b2
# ggsave("figure_5b-2.pdf", plot = p5b2, width=10, height=3, useDingbats=FALSE)
# ggsave("figure_5b-2.png", plot = p5b2, width=10, height=3, dpi=600)
#Plot major/minor allele ratio (Log2) with red/blue gradient
p5b3 <- ggplot(data=indic, aes(x=Manual_Rank_fig5, y=1, fill= Log2ratio)) +
geom_tile(color = "white")+
scale_fill_gradient2(low = "#4f94cd", high = "#ee2c2c", mid = "white",
midpoint = 0, limit = c(-2,2), space = "Lab") +
theme_minimal()+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, size = 12, hjust = 1),
axis.title.y=element_blank())+
coord_fixed()
p5b3
# ggsave("figure_5b-3.pdf", plot = p5b3, width=10, height=3, useDingbats=FALSE)
# ggsave("figure_5b-3.png", plot = p5b3, width=10, height=3, dpi=600)Supplemental Figure 1
#Supplemental Figure 1 color palette
farrowAndBall_palette <- c(
"#4d5b6a" #Stiffkey Blue
,"#6a90b4" #Cook's Blue
,"#599ec4" #ST Giles Blue
,"#a1c5c8" #Blue Ground
,"#7997a1" #Stone Blue
,"#427e83" #Vardo
,"#84b59c" #Arsenic
,"#919f70" #Yeabridge Green
,"#686a47" #Bancha
,"#c8bd83" #Churlish Green
,"#cb9e59" #India Yellow
,"#ecc363" #Babouche
,"#c57b67" #Red Earth
,"#d65f3d" #Charlotte's Locks
,"#a04344" #Incarnadine
,"#bf7a8f" #Rangwali
,"#8d8089" #Brassica
,"#50414c" #Pelt
,"#e5e0db" #Strong White
,"#444546") #Off-BlackSupplemental Figure 1 (200 lines)
#import data from current working directory
rar_current <- readRDS("fig3_200line.rds") #same RDS that was used in figure3
rar_current## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 1189 taxa and 598 samples ]
## sample_data() Sample Data: [ 598 samples by 16 sample variables ]
## tax_table() Taxonomy Table: [ 1189 taxa by 7 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 1189 tips and 1188 internal nodes ]rar_trans <- transform_sample_counts(rar_current, function(x) 100 * x / sum(x) )
# Agglomerate Taxa
glom_Order <- tax_glom(rar_trans, taxrank = "Order")
glom_Order## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 100 taxa and 598 samples ]
## sample_data() Sample Data: [ 598 samples by 16 sample variables ]
## tax_table() Taxonomy Table: [ 100 taxa by 7 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 100 tips and 99 internal nodes ]## pick a taxonomic level
taxrank_level = "Order"
gloms <- get(paste0("glom_", taxrank_level))
dim(otu_table(gloms))[1]## [1] 100gloms## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 100 taxa and 598 samples ]
## sample_data() Sample Data: [ 598 samples by 16 sample variables ]
## tax_table() Taxonomy Table: [ 100 taxa by 7 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 100 tips and 99 internal nodes ]## get top 15 orders
others <- names(sort(taxa_sums(gloms), decreasing = TRUE)[16:length(names(taxa_sums(gloms)))])
gloms <- merge_taxa(gloms, others, 1)
gloms## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 16 taxa and 598 samples ]
## sample_data() Sample Data: [ 598 samples by 16 sample variables ]
## tax_table() Taxonomy Table: [ 16 taxa by 7 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 16 tips and 15 internal nodes ]# sort orders based on relative abundance
# get otu table
ratab <- data.frame(otu_table(gloms))
dim(ratab)## [1] 16 598# calculate average relative abundance for all the samples
ratab <- data.frame(rowMeans(ratab))
colnames(ratab) <- "Relative_Abundance"
# merge otu table with taxonomy table
ratab <- merge(as.data.frame(gloms@tax_table@.Data), ratab, by = "row.names")
# sorting taxa by abundance from high to low for RA
sort0 <- c()
sort0 <- ratab[with(ratab, order(ratab$Relative_Abundance)),][[taxrank_level]]
sort0 <- as.vector(sort0)
sort0[is.na(sort0)] <- "Other"
sort0 <- sort0[! sort0 %in% "Other"]
sort0 <- c("Other", sort0)
sort0## [1] "Other" "Verrucomicrobiales" "Myxococcales"
## [4] "Solibacterales" "Xanthomonadales" "BetaproteobacteriaCL"
## [7] "Gemmatimonadales" "Sphingomonadales" "Acidobacteriales"
## [10] "Rhodospirillales" "AcidobacteriaPH" "Rhizobiales"
## [13] "Sphingobacteriales" "Bacillales" "Burkholderiales"
## [16] "Actinomycetales"## melt dataframe
melts <- psmelt(gloms)
# convert Class to a character vector from a factor
str(melts[[taxrank_level]])## Factor w/ 15 levels "Acidobacteriales",..: 3 NA NA NA NA NA NA NA 3 NA ...melts[[taxrank_level]] <- as.character(melts[[taxrank_level]])
# replace all NAs with Other
melts[[taxrank_level]][is.na(melts[[taxrank_level]])] <- "Other"
# convert Class back to a factor
melts[[taxrank_level]] <- as.factor(melts[[taxrank_level]])
levels(factor(melts[[taxrank_level]]))## [1] "Acidobacteriales" "AcidobacteriaPH" "Actinomycetales"
## [4] "Bacillales" "BetaproteobacteriaCL" "Burkholderiales"
## [7] "Gemmatimonadales" "Myxococcales" "Other"
## [10] "Rhizobiales" "Rhodospirillales" "Solibacterales"
## [13] "Sphingobacteriales" "Sphingomonadales" "Verrucomicrobiales"
## [16] "Xanthomonadales"melts$Order <- factor(melts$Order, levels=sort0)
# sort taxa
summary(melts)## OTU Sample Abundance SampleID
## Length:9568 Length:9568 Min. : 0.4778 B4_sb100_H: 16
## Class :character Class :character 1st Qu.: 3.5722 B4_sb101_H: 16
## Mode :character Mode :character Median : 4.9667 B4_sb106_H: 16
## Mean : 6.2500 B4_sb108_H: 16
## 3rd Qu.: 6.5111 B4_sb109_H: 16
## Max. :41.2500 B4_sb11_H : 16
## (Other) :9472
## Modified.ID x.cor y.cor Column
## B4.sb100.H: 16 Min. : 0.50 Min. :-198.50 C6 :1600
## B4.sb101.H: 16 1st Qu.: 9.50 1st Qu.:-124.50 C7 :1600
## B4.sb106.H: 16 Median :20.50 Median : -74.50 C8 :1600
## B4.sb108.H: 16 Mean :17.68 Mean : -82.73 C5 :1584
## B4.sb109.H: 16 3rd Qu.:25.50 3rd Qu.: -36.50 C1 : 800
## B4.sb11.H : 16 Max. :29.50 Max. : -0.50 C2 : 800
## (Other) :9472 (Other):1584
## Row Block Line PI.number Prep.time Seq.time
## R1 : 256 B4:3184 sb100 : 48 PI533839: 96 early:1136 L1:1136
## R10 : 256 B5:3200 sb101 : 48 PI48770 : 48 late :8432 L2:2800
## R11 : 256 B6:3184 sb106 : 48 PI533750: 48 L3:5632
## R12 : 256 sb108 : 48 PI533776: 48
## R13 : 256 sb109 : 48 PI533789: 48
## R14 : 256 sb11 : 48 PI533800: 48
## (Other):8032 (Other):9280 (Other) :9232
## Harvest Sampletype Redo Group Rep
## 9/22/16:3184 Rhizosphere:9568 No:9568 :9568 Min. : NA
## 9/23/16:6384 1st Qu.: NA
## Median : NA
## Mean :NaN
## 3rd Qu.: NA
## Max. : NA
## NA's :9568
## Kingdom Phylum Class
## Bacteria:9568 Proteobacteria:4186 Alphaproteobacteria:1794
## Acidobacteria :1794 AcidobacteriaPH :1196
## Actinobacteria: 598 Betaproteobacteria :1196
## Bacteroidetes : 598 ActinobacteriaPH : 598
## Firmicutes : 598 Bacilli : 598
## (Other) :1196 (Other) :3588
## NA's : 598 NA's : 598
## Order
## Other : 598
## Verrucomicrobiales : 598
## Myxococcales : 598
## Solibacterales : 598
## Xanthomonadales : 598
## BetaproteobacteriaCL: 598
## (Other) :5980Plot 200 lines
# plot relative abundance for rhizosphere (200 lines, merged)
ps1_1 <- ggplot(data=melts, aes(x=Sampletype, y=Abundance, fill=factor(get(taxrank_level)))) +
geom_bar(stat="identity",position = "stack") +
scale_fill_manual(values = farrowAndBall_palette[1:16]) +
theme(axis.text.x=element_text(size=12,color="black",angle=90),
axis.text.y=element_text(size=12,color="black"),
axis.title=element_text(size=16,face="bold"),
text=element_text(size=16)) +
guides(fill=guide_legend(title=taxrank_level))
ps1_1
# ggsave("figure_S1-1.pdf", plot = ps1_1, width=4.5, height=8, useDingbats=FALSE)
# ggsave("figure_S1-1.png", plot = ps1_1, width=4.5, height=8, dpi=600)
# plot relative abundance for rhizosphere (200 lines, separate)
ps1_2 <-ggplot(data=melts, aes(x=Line, y=Abundance, fill=factor(get(taxrank_level)))) +
geom_bar(stat="identity", width=1, position = "fill") +
scale_fill_manual(values = farrowAndBall_palette[1:16]) +
theme(axis.text.x=element_text(size=12,color="black",angle=90),
axis.text.y=element_text(size=12,color="black"),
axis.title=element_text(size=16,face="bold"),
text=element_text(size=16)) +
guides(fill=guide_legend(title=taxrank_level))
ps1_2
# ggsave("figure_S1-2.pdf", plot = ps1_2, width=8, height=8, useDingbats=FALSE)
# ggsave("figure_S1-2.png", plot = ps1_2, width=8, height=8, dpi=600)Supplemental Figure 1 (24 lines) Note:24 line chunk must be run after the 200 line chunk above, due to dependencies
#import data from current working directory
rar_current <- readRDS("fig1_24line.rds") #same RDS that was used in figure1
rar_current ## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 1787 taxa and 204 samples ]
## sample_data() Sample Data: [ 204 samples by 5 sample variables ]
## tax_table() Taxonomy Table: [ 1787 taxa by 6 taxonomic ranks ]rar_trans <- transform_sample_counts(rar_current, function(x) 100 * x / sum(x) )
#Agglomerate Taxa
## pick a taxanomic level
taxrank_level = "Order"
gloms <- tax_glom(rar_trans, taxrank = taxrank_level)
gloms## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 101 taxa and 204 samples ]
## sample_data() Sample Data: [ 204 samples by 5 sample variables ]
## tax_table() Taxonomy Table: [ 101 taxa by 6 taxonomic ranks ]## get tax table
tax_24 <- as.data.frame(gloms@tax_table@.Data)
others <- rownames(subset(tax_24, ! tax_24$Order %in% sort0))
## get 15 orders
gloms <- merge_taxa(gloms, others, 1)
gloms## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 16 taxa and 204 samples ]
## sample_data() Sample Data: [ 204 samples by 5 sample variables ]
## tax_table() Taxonomy Table: [ 16 taxa by 6 taxonomic ranks ]## melt dataframe
melts <- psmelt(gloms)
#convert Class to a character vector from a factor
str(melts[[taxrank_level]])## Factor w/ 15 levels "Acidobacteriales",..: 3 3 3 3 13 3 13 NA 3 3 ...melts[[taxrank_level]] <- as.character(melts[[taxrank_level]])
#replace all NAs with Other
melts[[taxrank_level]][is.na(melts[[taxrank_level]])] <- "Other"
# convert Class back to a factor
melts[[taxrank_level]] <- as.factor(melts[[taxrank_level]])
levels(factor(melts[[taxrank_level]]))## [1] "Acidobacteriales" "AcidobacteriaPH" "Actinomycetales"
## [4] "Bacillales" "BetaproteobacteriaCL" "Burkholderiales"
## [7] "Gemmatimonadales" "Myxococcales" "Other"
## [10] "Rhizobiales" "Rhodospirillales" "Solibacterales"
## [13] "Sphingobacteriales" "Sphingomonadales" "Verrucomicrobiales"
## [16] "Xanthomonadales"melts$Order <- factor(melts$Order, levels=sort0)
#sort taxa
summary(melts)## OTU Sample Abundance SampleID
## Length:3264 Length:3264 Min. : 0.000 B4_112_Leaf : 16
## Class :character Class :character 1st Qu.: 0.240 B4_112_Rhizo: 16
## Mode :character Mode :character Median : 2.960 B4_112_Root : 16
## Mean : 4.412 B4_163_Leaf : 16
## 3rd Qu.: 5.640 B4_163_Rhizo: 16
## Max. :48.720 B4_163_Root : 16
## (Other) :3168
## Block Line Sampletype Replicate Kingdom
## B4:1152 sb112 : 144 Leaf : 960 Min. : NA Bacteria:3264
## B5:1072 sb163 : 144 Rhizo:1152 1st Qu.: NA
## B6:1040 sb196 : 144 Root :1152 Median : NA
## sb201 : 144 Mean :NaN
## sb258 : 144 3rd Qu.: NA
## sb316 : 144 Max. : NA
## (Other):2400 NA's :3264
## Phylum Class Order
## Proteobacteria:1428 Alphaproteobacteria: 612 Other : 204
## Acidobacteria : 612 AcidobacteriaPH : 408 Verrucomicrobiales : 204
## Actinobacteria: 204 Betaproteobacteria : 408 Myxococcales : 204
## Bacteroidetes : 204 ActinobacteriaPH : 204 Solibacterales : 204
## Firmicutes : 204 Bacilli : 204 Xanthomonadales : 204
## (Other) : 408 (Other) :1224 BetaproteobacteriaCL: 204
## NA's : 204 NA's : 204 (Other) :2040melts$Sampletype <- factor(melts$Sampletype, levels = c("Leaf", "Root", "Rhizo"))Plot 24 lines
ps1_3 <- ggplot(data=melts, aes(x=Sampletype, y=Abundance, fill=factor(get(taxrank_level)))) +
geom_bar(stat="identity",position = "fill") +
scale_fill_manual(values = farrowAndBall_palette[1:16]) +
theme(axis.text.x=element_text(size=12,color="black",angle=90),
axis.text.y=element_text(size=12,color="black"),
axis.title=element_text(size=16,face="bold"),
text=element_text(size=16)) +
guides(fill=guide_legend(title=taxrank_level))
ps1_3
# ggsave("figure_S1-3.pdf", plot = ps1_3, width=6, height=8, useDingbats=FALSE)
# ggsave("figure_S1-3.png", plot = ps1_3, width=6, height=8, dpi=600)Supplemental Figure2
Supplemental Figure 2a
#import data from current working directory
pcs <- read.table("figS2_top10PC.txt", header = TRUE)
pcs$Top10PCs <- factor(pcs$Top10PCs, levels = c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6", "PC7", "PC8", "PC9", "PC10"))Plot for Supplemental Figure 2a
ps2a <- ggplot(data=pcs, aes(x=Top10PCs, y=ProportionExplained, fill=level)) +
geom_bar(stat="identity") +
scale_fill_manual(values=c("#E09E19", "#2D637F"))
ps2a
# ggsave("figure_S2a.pdf", plot = ps2a, width=6, height=4, useDingbats=FALSE)
# ggsave("figure_S2a.png", plot = ps2a, width=6, height=4, dpi=600)Supplemental Figure 2b
#get snp region of interest for chromosome4 (used to specify SNPs within ch4 locus)
qval <- read.csv("PC1.assoc.txt", header = TRUE, sep = "\t")
ch4 <- subset(qval, chr == 4)
goi <- subset(ch4, ps > 47500000 & ps < 49000000)
snpsOfInterest4 <- levels(factor(goi$rs))
#get snp of interest for chromosome6 (used to specify SNPs within ch6 locus)
ch6 <- subset(qval, chr == 6)
goi <- subset(ch6, ps > 50000000 & ps < 51000000)
snpsOfInterest6 <- levels(factor(goi$rs))subset for PC1 / chromosome 4
qval <- read.csv("PC1.assoc.txt", header = TRUE, sep = "\t")
qval$log10p <- -log10(qval$p_score)
#qval$p_wald_fdr <- p.adjust(a$p_wald, method="fdr") ##confirm for multi test correnction
chr4 <- subset(qval, chr ==4)
qval <- chr4
geno <- data.frame(CHR=qval$chr, BP=qval$ps, P=qval$p_score, SNP=qval$rs)Plot for Supplemental Figure 2b (PC1)
ps2b1 <- qqman::manhattan(na.omit(geno), ylim=c(0,4.5), suggestiveline = F, col = c("#2D637F"), highlight = snpsOfInterest4) #export pdf size width=3, height=6, png width = 300, height = 600
#qqfun(na.omit(geno), ylim=c(0,4.5), suggestiveline = F, col = c("#2D637F"), highlight = snpsOfInterest4) #NOTE: "qqfun" = "qqman::manhattan", with source code to change SNP color from green to red.subset for PC5 / chromosome 6
qval <- read.csv("PC5.assoc.txt", header = TRUE, sep = "\t")
qval$log10p <- -log10(qval$p_score)
#qval$p_wald_fdr <- p.adjust(a$p_wald, method="fdr") ##confrim for multi test correnction
chr6 <- subset(qval, chr ==6)
qval <- chr6
geno <- data.frame(CHR=qval$chr, BP=qval$ps, P=qval$p_score, SNP=qval$rs)Plot for Supplemental Figure 2b (PC5)
ps2b2 <- qqman::manhattan(na.omit(geno), ylim=c(0,4.5), suggestiveline = F, col = c("#2D637F"), highlight = snpsOfInterest6) #export pdf size width=3, height=6, png width = 300, height = 600
#qqfun(na.omit(geno), ylim=c(0,4.5), suggestiveline = F, col = c("#2D637F"), highlight = snpsOfInterest6) #NOTE: "qqfun" = "qqman::manhattan", with source code to change SNP color from green to red. subset for PC10 / chromosome 6
qval <- read.csv("PC10.assoc.txt", header = TRUE, sep = "\t")
qval$log10p <- -log10(qval$p_score)
#qval$p_wald_fdr <- p.adjust(a$p_wald, method="fdr") ##confrim for multi test correnction
chr6 <- subset(qval, chr ==6)
qval <- chr6
geno <- data.frame(CHR=qval$chr, BP=qval$ps, P=qval$p_score, SNP=qval$rs)Plot for Supplemental Figure 2b (PC10)
ps2b3 <- qqman::manhattan(na.omit(geno), ylim=c(0,4.5), suggestiveline = F, col = c("#2D637F"), highlight = snpsOfInterest6) #export pdf size width=3, height=6, png width = 300, height = 600
#qqfun(na.omit(geno), ylim=c(0,4.5), suggestiveline = F, col = c("#2D637F"), highlight = snpsOfInterest6) #NOTE: "qqfun" = "qqman::manhattan", with source code to change SNP color from green to red. Supplemental figure 2c
#OTU list generated from Figure4_get_otus.R
#the list is called out_all
#OTU within the chromosome4 locus
out0 <- rev(c("PC1","X692","X14","X172","X544","X798","X1170","X78","X30","X312","X766","X1020","X199","X734","X318","X577",
"X947","X1190","X190","X445","X302","X531","X1289","X11","X325","X1239","X904","X943","X1039","X506","X94",
"X82","X326","X461","X695","X286","X205","X148","X340","X996"))
#OTU within for chromosome6 locus
out1 <- c("PC10","PC5","X323","X1103","X802","X359","X846","X422","X730","X521","X682","X1158")
#OTU shared with both chromosome4 and 6 loci
out_all <- c(out1[! out1 %in% "X18"], "X18", out0[! out0 %in% "X18"])
otuhits <- out_allSupplemental figure 2c (Ch4)
#get snps data
goi <- data.frame()
for(i in 1:length(otuhits)){
oid <- otuhits[i]
print(i)
print(oid)
file <- paste0(oid, ".assoc.txt")
a <- fread(file, data.table=FALSE)
ch4 <- subset(a, chr == 4)
ch4$file <- file
ch4$otu <- gsub(".*\\/|.assoc.txt", "", ch4$file)
goi_add <- subset(ch4, ps > 47500000 & ps < 49000000)
goi <- rbind(goi, goi_add)
}
#rank otu
goi$otu <- factor(x = goi$otu, levels = out_all, ordered = TRUE)
length(levels(factor(goi$rs)))Plot for Supplemental Figure 2c (ch4)
ps2c1 <- ggplot(data = goi, aes(x = rs, y = otu)) +
geom_tile(aes(fill = -log10(goi$p_score))) +
scale_fill_gradientn(colours=c("white","#e5e0db","#6a90b4","#2D637F","#4d5b6a"), values = rescale(c(0,1.5,3,6,9)),guide = "colorbar") +
theme(axis.title.y = element_blank(),
axis.ticks.y = element_blank(),
legend.position = "top")
ps2c1
# ggsave("figure_S2c-1.pdf", plot = ps2c1, width=6, height=8, useDingbats=FALSE)
# ggsave("figure_S2c-1.png", plot = ps2c1, width=6, height=8, dpi=600)Supplemental figure 2c (Ch6)
#get snps data
goi <- data.frame()
for(i in 1:length(otuhits)){
oid <- otuhits[i]
print(i)
print(oid)
file <- paste0(oid, ".assoc.txt")
a <- fread(file, data.table=FALSE)
ch6 <- subset(a, chr == 6)
ch6$file <- file
ch6$otu <- gsub(".*\\/|.assoc.txt", "", ch6$file)
goi_add <- subset(ch6, ps > 50000000 & ps < 51000000)
goi <- rbind(goi, goi_add)
}
##rank otu
goi$otu <- factor(x = goi$otu, levels = out_all, ordered = TRUE)
length(levels(factor(goi$rs)))Plot for Supplemental Figure 2c (ch6)
ps2c2 <- ggplot(data = goi, aes(x = rs, y = otu)) +
geom_tile(aes(fill = -log10(goi$p_score))) +
scale_fill_gradientn(colours=c("white","#e5e0db","#6a90b4","#2D637F","#4d5b6a"), values = rescale(c(0,1.5,3,6,9)),guide = "colorbar") +
theme(axis.title.y = element_blank(),
axis.ticks.y = element_blank(),
legend.position = "top")
# ggsave("figure_S2c-2.pdf", plot = ps2c2, width=6, height=8, useDingbats=FALSE)
# ggsave("figure_S2c-2.png", plot = ps2c2, width=6, height=8, dpi=600)Supplemental figure 2c (bacterial order)
#import data from current working directory
qval <-read.csv("figS2_order.txt")
#Color palette used for supplemental figure 2c
FS2C_palette <- c("#cb9e59"#"Acidobacteriales"
,"#a1c5c8"#"Actinomycetales"
,"#8B0000"#"Bacillales"
,"#c8bd83"#"Burkholderiales"
,"#84b59c"#"Gemmatimonadales"
,"#427e83"#"Myxococcales"
,"#444546" #"Other"
,"#ecc363" #Rhizobiales
,"#919f70"#Rhodospirillales
,"#a04344"#"Solibacterales"
,"#c57b67"#"Solirubrobacterales"
,"#bf7a8f"#"Spartobacteriales"
,"#599ec4"#"Sphingobacteriales"
,"#7997a1"#"Verrucomicrobiales"
,"#d65f3d")#"Xanthomonadales" Plot for Supplemental Figure 2c (bacterial order)
ps2c3 <- ggplot(data=qval, aes(x=rank, y=1, fill= Order)) +
geom_tile()+
theme_minimal()+
scale_fill_manual(values = (FS2C_palette)) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, size = 12, hjust = 1),
axis.title.y=element_blank())+
coord_fixed()
ps2c3
# ggsave("figure_S2c-3.pdf", plot = ps2c3, width=10, height=6, useDingbats=FALSE)
# ggsave("figure_S2c-3.png", plot = ps2c3, width=10, height=6, dpi=600)